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ABSTRACT 



An algorithm for the computer-aided design of passive electrical 
filters using minimization techniques is presented. The configuration 
or topology of the filter and its desired frequency response are speci 
fied. A measure of the deviation between the desired and the actual 
frequency response is minimized using the gradient projection search 
method. Examples are presented which demonstrate the performance of 
the minimization procedure. 
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Io INTRODUCTION 



The development of extremely high-speed computers and efficient 
algorithms has opened the door to the use of iteration techniques that 
heretofore have been prohibitively costly in time and/or money. This 
high-speed capability makes it possible to consider a computer-aided 
solution to the problem of designing an electrical filter to match — 
as closely as possible — a given frequency responses wherein many 
hundreds — sometimes thousands — of successive solutions of the network 
frequency response are required. Use of the designing engineer's expe- 
rience in adaptation of filters whose response is already known, but is 
not quite that desired, usually demands laborious hand calculation. In 
addition, there is occasional difficulty in achieving the required 
"closeness of fit" to the desired response. Some frequency response 
requirements may lie beyond the capability of synthesis techniques, 
being solvable only by experience coupled with trial and error proce- 
dures — which may never produce the desired accuracy. 

With the availability of a computer program capable of systematically 
minimizing deviations from the required response, design time can be cut 
dramatically. The ease with which this tool can produce results has 
multiple benefits: 

1. The engineer is free to concentrate on design without the tedium 
of repetitious computation. 

2. The designer can investigate numerous circuit configurations to 
determine one which "best" fits the needs of simplicity, cost, and 
quality of performance. 

3. The computer's speed enables the designer to quickly obtain an 
answer. 
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II. OPTIMIZATION THEORY 



A, MINIMIZING THROUGH OPTIMIZATION 

Optimization techniques are used to provide the means of systemat- 
ically and effectively developing the solution to the problem. Classical 
methods for designing filters to meet a given frequency response are 
highly time-consuming, and usually produce only approximations of the 
desired response, perhaps matching the general shape, and a specified 
bandwidth [Ref, 1], Optimization techniques utilize the designer's 
particular definition for "closeness of fit" of the designed to the 
specified frequency response and then proceed through successive adjust- 
ments of the variable parameters to minimize the deviation. The function 
defined by the designer to specify this closeness of fit will be referred 
to as "the cost function", 

A specific filter configuration is required by the optimization 
process; no provision is made to permit computer-determined alteration of 
the circuit in an effort to improve the overall closeness of fit. This 
decision is left to the engineer to be carried out during his interaction 
with the program, 

B, MINIMIZATION TECHNIQUES 

Three methods for minimization were considered, direct, or pattern 
search; the gradient method; and second-order methods, 

1 , Pattern Searc h 

The pattern-search method makes use of successive perturbations 
of the variable parameters to determine when a lower value of the cost 
function is achieved. When a reduction in value of the cost function is 
determined, the new point is then used to again perturb the parameters, 
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taking larger steps when continued successes are achieved. If a failure 
to decrease the cost function occurs, the step size is reduced until 
success is achieved. This is the simplest of the three methods, re- 
quiring only values of the cost function, but is somewhat sensitive to 
starting point when local minima are present [Ref, 2], 

2 , Gradient Meth od 

The gradient method makes use of the negative of the first 
partial derivative of the cost function with respect to each of the 
variable parameters; this determines the direction of steepest descent, 

A step of specified distance is taken in the steepest-descent direction 
to decrease the value of the cost function, 

A closely related method, that of gradient projection [Ref, 3], 
makes use of the negative gradient at the starting point to extend a 
vector in this "steepest-descent" direction until one of the parameter 
constraints is reached. The smallest cost function value along this 
extended vector is then determined. The negative gradient vector is 
computed for this new point, and the previous step is repeated until the 
minimum is reached. The algorithm will always move in the negative 
gradient direction unless doing so causes a violation of the constraints, 
in which case the negative gradient vector is projected onto the con- 
straint boundary, and the next step is taken in the direction of this 
projection, 

3. Second-Order Methods 

Second-order methods involve the use of both the first and 
second partial derivatives of the cost function, and were not considered 
once the difficulty in generating the second partial derivative for 
filter networks was recognized 
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C. PROBLEM CHARACTERISTICS 
1. The Cost Function 



The general format for the cost function, J(p), used in the 
optimization scheme is defined as: 

0 ) 

max N 

J(p) * / h(m) |e(p)| N d u (2.1) 

“min 

Here p is the vector of variable circuit parameters; h (to) is a weighting 
function used to force the solution to more closely fit the prescribed 
response in certain critical areas; e(p) is the error between the desired 
and the actual frequency response; N represents the power to which e(p) 
is raised to determine the type of closeness of fit (e.g., N = 2 for mean- 

square), and normally takes on the value of 1 , 2 or 4; a> m ‘ n and u are 

the lower and upper frequencies over which the integration is performed. 

The integration must be written in discrete form for use in a 
computer program; thus the following approximation is used; 

K 

J(p) = Ad) £ h(ti) • +kAoo) i R,l(<»> .s + kAll))-R a ( a, min + ^Au) | ^ . (2.2) 

V. min 1 d min a min 1 

In this equation Au> is the difference between successive frequency points; 

K represents the number of frequency points to be matched; N is as 

described above, and h(u> . +kAw) is the discrete equivalent of h(w); 

min 

t h 

R d^“min + * <A “^ 1S desired frequency response at the k frequency 

point, while R (oj . +kAu) is the actual response. It is this function, 
d mi n 

J(p), that the optimization process is to minimize. If a perfect fit to 
the desired frequency response is achieved, the value of J (p ) will be 
zero, its absolute minimum. 
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2. Termination Criteria 



Iteration of an optimization program can be terminated i i several 
ways. Three of these are: 

1. When further iteration produces no meaningful improvement in close- 
ness of fit to the desired frequency response, iteration is terminated. 

A suitable equation is: 



where e is some relatively small positive number. If e is selected to be 
too small, excessive iterations will usually be required to produce two 
successive values of J(p) close enough together to terminate the procedure, 
resulting in wasted computer time. If a large value is chosen for e , the 
program may terminate far from a minimum. 

2. When the norm of the gradient vector is quite close to zero, either 
a minimum, a maximum or a saddle point has been found; in a constrained 
problem, maximums are less likely. Numerous problems were run, and 
neither a maximum nor a saddle point ever resulted as an optimization 
solution. This criterion can be written: 



where y is a suitably small positive number. 

3. A third criterion, not written into the computer program used in 
this investigation, but one which might prove useful to the designer in 
utilization of the program is: 




(2.3) 



0(P) 



J^(p) < 6 



(2.5) 
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where 6 is some suitably small positive number. Inclusion of this 
inequality into the computer program would cause termination once the 
cost function reaches a value less than 6. 

3. Constraints 

In many practical design problems, there is some idea of a 
reasonable range of constraints for each variable parameter, therefore 
it is assumed that constraints of the form shown below can be imposed on 
the n variable parameters: 

a n 1 P-j 1 b -j i = l,2,...,n (2.6) 

where a = (a^ .a^, . . . ,a n ) represents the vector lower limit, and 
b = (b^ ,b 2 ». . . ,b ) is the vector upper limit. In the case of several 
variables, the constraints influence the rate of convergence, especially 
for the gradient-projection method. For example, a reduction of con- 
straints from 0.<p <1. to 0. < p < .5 , where n s 5, will reduce 
- K n - - 'n - 

C 

a five-dimensional space by 96.9 percent (using (1/2) - 03124), leaving 

a smaller region in which to search for the minimum 
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III. THE OPTIMIZATION PROGRAM 



A„ THE FREQUENCY ANALYSIS PROGRAM 

Several existing computer programs for calculating the frequency 

response of a given circuit were studied. Some of these were: 

Calahan 

Sceptre 

Ecap 

Cornap 

The Calahan Network Analysis Program was chosen as being best able to 
provide the required data, in the desired form, with high accuracy and 
reasonable speed, 

1 „ The Basic Program 

The Linear Network Analysis Program by D„ A, Calahan [Ref, 4] is 
a large general purpose program capable of producing a variety of circuit 
information. Inputs to the program are: either the network transfer 

function or a topological description; keys to indicate the desired output; 
high- and low-frequency limits and the number of points to be determined, 
if the frequency response is desired; initial and final values, initial 
and final time, and the time increment if a time response is desired. 
Provisions are made for varying one element for as many as eight differ- 
ent values. ' The output from the program includes the coefficients of the 
network function, the poles and zeros, the frequency response, the time 
response (to an arbitrary input), and the symbolic form of the network 
function. The program is capable of accepting circuits with as many as 
thirty nodes, one hundred passive elements, and twenty active elements. 
Printer-plot graphs of frequency and time response are produced when 
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desired. Negative element values, as well as negative frequencies, are 
acceptable to the program, 

2. Accuracy 

Examination of the accuracy of the frequency response produced 
by the Calahan analysis was accomplished in the following fashion. 
Specific circuits whose frequency response could oe readily calculated 
using formulae (Butterworth and Chebychev low-pass filters in particular 
were used, producing the expected frequency response using formulae 
found in [Ref, 5]) were tested using the nominal element values. The 
Calahan results compared favorably with the answers from the formulae. 
For example, analysis of a five-element Butterworth filter gave the 
following results: 





Formula 


Calahan 


Freq(Hz) 


Gain(db) 


Gain(db) 


,15 


-2,10196 


-2.10181 


.16 


-3.10327 


-3J0314 


.17 


-4 30471 


-4.30461 


.18 


-5,65472 


-5 65465 


,19 


-7,09729 


-7,09723 


.20 


-8,58441 


-8, 58439 


3. Modifications 


Because only the frequency response was 


desired as an output from 


Calahan, and this on a highly repetitive basis. 


modifications to the 


program were made. Several 


subroutines and the 


main program were removed; 



a new main program was substituted so that the solution to the minimization 
problem could be generated by repeated use of the Calahan subroutines to 
determine the required frequency response. Two of the remaining Calahan 
subroutines were modified to conform to the desired iteration scheme, 

B. THE MINIMIZATION PROGRAM 

The gradient, or steepest-descent method, was chosen to be used for 
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minimization of the previously specified cost function,, (The pattern 
search solution to this problem is the subject of a Naval Postgraduate 
School thesis by Captain James Lau, USMC. Comparative results are 
discussed in Chapter VI.) 

1 . Determination of the Gradient Vector 

Some thought was given to the best method of producing the vector 

90(p) 

— rr— . The literal matrix of partial derivatives was manually computed 
°P 

for three- and five-element filters, and it was quickly seen that this 
method is much too cumbersome and time-consuming, and will greatly reduce 
the ease with which the designer can use the program. As an example, 
given below is one partral for a three-element filter (two variable 
elements, L and C; R=l). 

3 J(L,C) = * -2[LC 2 (kA^) 4 +L(kA ( o) 2 -c(kAa ) ) 2 1 

9L k-1 [L 2 C 2 (kAijj) 4 +(L 2 -2LC) (kA(jj) 2 +l ] 2 

{___ _£ _ _ __ - — - — T-j » (3.1 ) 

L C (kAco) +(L -2LC) (kAw) +1 l + (kAio) 



It was then decided to investigate the method of finite-difference 
approximation. This method makes use of the definition of the partial 
derivative. 



aJ(p) 

ap-l 



1 im 
AP^O 



J^Pl+Ap-j ,P2> ° » » jP n ) ~ J ^ Pi j Pg » ° ° • ,P^j ) 
— 



(3.2) 



The finite-difference approximation assigns a small but finite value to 
AP-| and determines the approximate partial derivative as: 



9 J ( P ) J ( P"| ^"^P-| 9 P2 9 ° ° “ 9 P^ ) ~ J(P*j 9 P^ 9 " 0 ° 9 P^ ) 

__ . _____ -- 



(3.3) 
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Experimentation with various networks, using the Calahan subrou- 
tines to produce the frequency response indicated that a Ap of 0,01 
percent gives reasonably accurate results. This accuracy was determined 
in two ways: by comparison with the values of the partial derivative 

vector for a two-variable-element filter produced by literal equations, 
one of which is equation (3.1); by use of discrete approximation to 
produce the gradient vector at the known absolute minimum of various 
larger networks. The first method produced normalized discrete approx- 
imation vectors within 1.5 percent of the true partial derivative vector 
in each of the ten points compared. The second method utilized larger 
networks whose actual frequency response was known, as well as the element 
values necessary to produce this response. When the correct element 

values were entered into the optimization program, each element of the 

-4 

gradient vector was less than 5 x 10 , in all examples tested, 

2. The Gradient Method 

Initially, a number of circuit optimization problems were tried, 
using the gradient method. In this program, new values for the variables 
are produced using the following iteration equation: 



P 



(i+1) 




dJ(p ll) ) 

^ TTT " 



(3.4) 



Each iteration produces a new set of parameter values, where the normal- 
ized gradient vector is "stepped" a distance g. When the new value for 
the cost function is greater than the previous value, 6 is halved and the 
step is repeated until a smaller function value is obtained. This method 
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produces good results on filters with up to five variable parameters, but 
required many iterations that reduce the value of the cost function very 
slowly as the minimum is approached. A better method was therefore 
sought. 

2. The Gradient-Projection Method 

In hope of eliminating some of the disadvantages encountered in 
the gradient method, the gradient-projection method introduced by J. B. 
Rosen was investigated. This algorithm showed promise of more rapid 
convergence. Additionally, a program using this method was locally 
available. This is the Gradient Projection Program written by Dr. D. E. 
Kirk [Ref. 6]. The process is an iterative numerical method for locating 
the extremum of a function of many variables which are required to satisfy 
linear constraints. Unconstrained variables cannot be used with this 
method, as will be shown. 

A simple two-dimensional example will serve to illustrate the use 
of this method. Let f(p) be a differentiable function of two variables 
Pl and p^; f ( p^ , p^ ) represents the value of the function at the point 
(p-|,p 2 ). The objective is to find the point (P^*,P 2 *) where f(p*) is a 
minimum, and which satisfies the constraints: 

1 P] 1 

0. <_ p 2 <_ 7. 

Figure 1(a) represents a possible contour plot of f(p) as a function of 

p-| and p^; the region which satisfies the constraints lies within the 

rectangle bounded by the lines Cl, C2, C3 and C4. The starting point is 

at p(0) = (2,1). At this point, the negative gradient vector, - , 

9p 

is calculated, and extended until a constraint boundary is reached, in 
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Figure 1 

Contour Maps Demonstrating the Gradient Projection Method 
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this case C4 The gradient at p(l), the intersection of C4 and the ex- 



af(0) 



3f(0) 



tended vector — ^ , is calculated. The inner produce! of - 

0 P 0 P 

and " is negative, indicating that the minimum value of f(p) 

along this extended vector is somewhere between p(0) and p(l); therefore 
interpolation is necessary. Successive interpolation locates the point 
p(2), where - is orthogonal to the originally extended vector. 

Orthogonality is indicated when the inner product of - and - 

is zero. The negative gradient vector - ^ is then extended to the 

nearest constraint boundary, and the process continues until the minimum 
is reached at p(6), as shown in Figure 1(a). 

Figure 1(b) shows a different hypothetical contour plot for f(p), 
using the constraints: 



0. < p 1 <_ 10. 

0. <_ p^ 5. 



In this case, the absolute minimum lies outside the allowable constrained 
region; the minimum therefore lies along a boundary. This situation 
causes the constraint boundary C3 to be "active"; that is, it separates 
the acceptable region (that region within the constraint boundaries) from 
the absolute minimum of the function f(p). Once again, the starting point 
is p(0) = (2,1). As before, the inner product of - and - — ^ — 

is negative, indicating that interpolation is necessary. The point p(2) 
is found and its negative gradient extended to the boundary C3 where, 
unlike the preceeding step, the inner product of - anci ~ is 

positive, indicating that p(3) is the minimum along the extension of 
- and that no interpolation is necessary. However, a change of p 

in the direction of - would cause a violation of the constraints; 

dP 

therefore p is changed along the vector represented by the projection of 



19 



- along the boundary C3» This projection is denoted by P[- , 

and is extended until another constraint is intersected, this time C4. 

Here the inner product of P[- — and - is negative, requiring 

interpolation,. Successive interpolations locate p(5), where P[- 

ap 

and - are orthogonal. This is the minimum within the acceptable 

ap 

region. It should now be evident that unconstrained parameters cannot 
be used with this method. 

In Figures 1(a) and (b), the functions to be minimized are uni- 
modal. If, however, a given function is not unimodal within the acceptable 
region, then more than one minimum can exist, and the minimum found by the 
gradient-projection method may depend on the starting point. Figure 2 
demonstrates this, showing solutions for two different starting points 
(gradient vectors are not labeled in order to enhance clarity). A point 
of interest is demonstrated by the projection vector from the second 
starting point: the actual minimum along this projection vector is 

located at A; however, the inner-product method indicates that no inter- 
polation is necessary, and accepts point p'(l) as the next starting point. 
This clearly demonstrates the requirement of unimodal ity for a function 
to be minimized by this interpolation method if the absolute minimum is 
desired, and multiple starting points are to be avoided. 

A more detailed explanation of this method, including the matrix 
equations required to be solved in the computer program, is contained in 
Reference 6. 

4. Cost-Function Complexity 

The inner-product interpolation scheme incorporated in the gradient 
projection program produces accurate results when a convex function is To 
be minimized; a convex function is unimodal. In using this interpolation 
method to optimize the frequency response of a filter network, it was found 
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Contour Map Showing Location of Two Different 
Minima by the Gradient Projection Method 



that the scheme was not functioning so that even-decreasing values of the 
cost function were obtained. This helped confirm what was already sus- 
pected: the cost function, J(p), as applied to filter network frequency 

response may not be a convex function. An actual contour map of 0(p) as 
a function of two variable parameters is shown in Figures 3(a) and (b). 

This map was produced oy solving for the cost function over a grid of 
values of the two variables and joining points of equal cost-function 
value. The map shows that while only one minimum exists, the function 
0 ( p ) is obviously not convex. Higher-dimensional maps cannot be easily 
drawn, but experience in minimizing filters with several variable elements 
indicates that more than one minimum may be present. Therefore it is 
concluded that the cost function, J(p), as used in this procedure is not 
convex, and is possibly, in general, not ummodal. No conclusive results 
in this area were sought, as this was not the object of the investigation. 

5. The Function-Search Interpolation Procedure 

Armed with the knowledge of the cost function's complexity, a new 
interpolation method was sought. This method would ideally achieve the 
following goals: 

Goal 1 Continually decrease the value of the cost function, 
never accepting any new point unless it yields a lower value of J(p). 

Goal 2, Make efficient use of computer time. 

Goal 3, Increase the likelihood of converging to the absolute 
minimum, if local minima should be present 

Investigation of various possible methods produced one that showed 
excellent promise of fulfilling all three requirements. This interpolation 
does not use successive values of the gradient vector, but successive 
values of the cost function A simple example is the easiest means of 
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:r=1j2 




Figure 3 

Actual Contour Map Showing Lines of Equal Value of 
J(p) as a Function of the Two Variable Elements 
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22 




C=p 



1 

1 



(a) 



!R=ln 




Figure 3 

Actual Contour Map Showing Lines of Equal Value of 
J ( p) as a Function of the Two Variable Elements 
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explaining the interpolation procedure. Using a filter with two variable 
elements (so that geometrical representations are possible), a contour 
map of J(p) as a function of the variables might be as shown in Figure 4. 




Hypothetical Contour Map Showing the Segmented 
Extension of the Negative Gradient Vector 

At the base point where the cost is 0(1), the negative gradient is calcu- 
lated and extended until it meets the first constraint boundary, p-|(min), 
as would have been done with the inner product interpolation method. The 
length of the extended vector will be referred to as T . Now, however, 
the line in the direction of the extended vector is subdivided into an 
even number of intervals (four intervals and five points will be used in 
this example). The cost function is determined for each of the remaining 
points, J(2) through J(5). The smallest value of J and its corresponding 
point number are determined. If J(5) is the smallest cost-function value, 
no interpolation is performed; if J(l) is the smallest, the interval 
between J(l) and J(2) is investigated, since the gradient vector's direc- 
tion indicates that some point with a smaller value than J(l) lies in this 
direction. If an interior point is smallest, interpolation must be * 
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accomplished. As is seen in Figure 4, 0(2) has the smallest cost-function 
value. The interpolation scheme then passes a parabola through 0(2) and 
its adjacent points (J(l) and J(3)), determines the minimum point of this 
parabola along the line in the direction of the extended vector, T, and 
evaluates the cost function, J(6), at this new point, as shown in Figure 
5(a). The purpose of approximating the actual curve of the cost function 
by a parabola, is that for reasonably well-behaved functions, a parabola 
joining three relatively close points on the curve will closely approxi- 
mate the actual curve. This helps to increase the speed with which the 
actual minimum is found. 

A detailed explanation of the procedure for finding the parabolic 
minimum is contained in Appendix A. 

Having solved for J(6) as shown in Figure 5(a), the points J(2), 

J(6), and J(3) are then used to solve for another parabolic minimum 

(Figure 5(b)). This process continues until the difference in extended 

vector length of the location of two successive minima is within five 

percent of T . In this case, location of parabolic minima would termi- 
max 

nate if T-, - T, = 0.05 T . The interpolation is then considered 
/ 6 max 

complete, and the lowest cost function value achieved is used as the new 
starting point. 

In numerous examples, this interpolation process has succeeded in 
quite reasonably accomplishing each of the desired goals. 

Goal 1. The requirement to continually decrease the value of the 
cost function was achieved. 

Goal 2. Computation time has been reduced over the inner* product 
interpolation method. Approximately the same number of points along the 
gradient vector are necessary in the function -search method, but only the 
value of the cost function, and not the gradient vector, is required at 
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Actual curve 
Parabolic estimate 
Actual minimum 
Parabolic minimum 



Distance Along the 
Extended Vector, T 




Distance Along the 
Extended Vector, T 



(b) 



Figure 5 

Parabolic Estimation Diagrams 
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each point. This eliminates n(m - 2) cost- function calculations, where 
n is the number of variables, and m is the number of points considered 
along the projection vector (the function search method still requires 
the gradient vector at the first and last points, accounting for the -2). 

Goal 3. The likelihood of converging to the absolute minimum 
appears to have been increased by use of this method. Results from its 
use in circuit frequency-response minimization indicate that such may be 
the case. Appendix B cites two examples (many more have been noted) 
where a distinct "break" occurred in what appeared to be the final approach 
stages where the gradient procedure converges to a minimum point. A sudden 
and reasonably large change in parameter values and decrease in cost- 
function value take place, something that does not occur using gradient 
methods that do not make use of the extended vector. 
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IV, RESULTS 



A general result is first presented, followed by specific examples 
illustrating results for several of the design problems investigated. 

The final section of the chapter discusses a problem encountered during 
the investigation, 

A. GENERAL RESULTS 

Since only a limited number of filters were investigated, the results 
are not conclusive but serve to demonstrate that optimization of the 
frequency response of a given filter is feasible using the techniques 
previously described. Although all the examples are of the passive variety, 
use of active elements should be feasible using models for the active 
elements required for the Calahan subroutines described in Reference 3. 

B. SPECIFIC RESULTS 

In this section, specific results for those filters investigated that 
representatively demonstrate the capability of the optimization program 
are presented with circuit diagrams, program solutions and graphical 
results. 

1 . Example 1 

Figure 6 contains the results of the optimization of a six-element 
low-pass filter. The desired curve is the Butterworth filter response 
[Ref. 1]: 

| Gain | = 20 Log 1Q [(1 + u> 2n )" 1/2 ] (4.1) 

where n is the number of reactive elements. In this case, n = 5. The 
actual response deviates from the desired only slightly. A constant 
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Circuit Configuration 




Starting Values of the Variable Parameters 



L ] = 1.0 L 2 = 1.1 L 3 = 1.0 C 1 = 2.4 C 2 = 0.9 



Final Values of the Variable Parameters 
L 1 = 1.557 L 2 = 1.695 l_ 3 = .500 C ] = 1.478 C 2 = 1.064 

Number of Iterations — 43 Computer Run-time - 136 seconds 

Termination Caused by - J^ 1+ ^ < 10 ^ 

Graphical Comparison of Desired and Actual Frequency Response 




f (Hz) 



Figure 6 

Data for Example 1 
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weighting matrix (h. = 1., i - 1,...,K) was used in this example, and in 
all other examples unless otherwise stated. 

Figure 7 demonstrates the rapidity with which the program improves 
the actual response. The filter response is shown in comparison to the 
desired response after the indicated number of iterations for the problem 
described in Figure 6. After ten iterations the actual response is a 
relatively good approximation to the desired response. 

2. Example 2 

Figure 8 contains the results of the optimization of a five- 
element low-pass filter, using straight-line approximation for the desired 
frequency response. Since this is not a realizeable response, it is only 
possible to approximate it with the actual response. The results are 
reasonably close to the desired shape and bandwidth. If it were desired 
to more closely match the 3-db bandwidth, a weighting matrix could be 
used that places greater emphasis on the points where the desired response 
is -3 decibels. 

3. Example 3 

Figure 9 contains the results of the optimization of an eight- 
element band-pass filter, wherein six elements were varied. The desired 
response was realizeable with the specified network configuration, and 
was very closely approximated, except in the frequency interval (0.14,0.22) 
where the deviation was noticeable. The bandwidth and general shape are 
reasonably accurate. This same filter circuit was optimized using a 
weighting matrix where both the upper and lower 3-db points were weighted 
with a value five times that of the weight of the other points. The devi- 
ation at the 3-db points was improved by a factor of about three. 
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Figure 7 

Graphs Showing Sequential Improvement in Actual Response 



Circuit Configuration 




Starting Values of the Variable Parameters 



2.0 



2.0 



C 1 = 2.0 



C 2 = 2.0 



’I 2 

Final Values of the Variable Parameters 



L 1 = 1.434 



L 2 = 1.093 



Number of Iterations - 40 



C 1 = 1.635 



0.385 



Computer Run-time — 132 seconds 



Termination Caused by j(i) _ j(i+l) < i 0 ”^ 



Graphical Comparison of Desired and Actual Frequency Response 




f- (Hz) 



Data for Example 2 
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Circuit Configuration 




Starting Values of the Variable Parameters 
L 1 = 1.0 L = 1.0 L 3 = 1.0 C-j =1.2 C 2 = 1.0 c 3 = 1.2 

Final Values of the Variable Parameters 
L ] = 0.231 L 2 = 0.520 L 3 = 0.237 C ] = 1.140 C 2 = 0.485 C 3 = 1.140 

Number of Iterations — 58 Computer Run-time — 328 seconds 

Termination Caused by - J^ 1+ ^ < 10~® 

Graphical Comparison of Desired and Actual Frequency Response 




Figure 9 

Data for Example 3 
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4. Example 4 



Figure 10 contains the results of the optimization of the same 
circuit used in example 3, this time specifying a straight-line approxi- 
mation to the desired frequency response,, Again, an approximate solution 
is the best that can be achieved, The final curve matches the desired 
reasonably well in general shape and bandwidth. 

5. Example 5 

Figure 11 contains the results of the optimization of a five- 
element band-pass filter where non-ideal inductors are considered. The 
desired frequency response is realizeable with the given circuit config- 
uration using ideal elements. When a resistance of R $ = 0.1 is included 
in series with both inductors, the optimized response is slightly lower 
in gain than desired, but still is a good approximation as shown in the 
graphical comparison. The same circuit was optimized with R g = 0.01 and 
with R s = 0.5. The former gave virtually an exact match to the desired 
response; the latter produced a response that was generally 3.5 decibels 
below the flat portion of the desired response, and somewhat degraded in 
shape. 

6. Example 6 

Figure 12 contains the results of the optimization of an eleven- 
element low-pass filter whose desired frequency response is described by 
equation (4.1), with n = 10, The optimization produced a close fit to 
the desired response (maximum deviation from the desired response was 0.27 
decibels), but at the expense of relatively long computer run-time. This 
example required slightly over 32 minutes and 99 iterations to produce 
the final solution. After 60 iterations, the same general response had 
been found, with a maximum deviation of 0.42decibels (and an approximate 
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Circuit Configuration 



C 3 L 3 




Starting Values of the Variable Parameters 
L-j = 1.0 L 2 = 1.0 L 3 = 1.0 C ] = 1.2 C 2 = 1.0 C 3 = 1.2 
Final Values of the Variable Parameters 
L ] = 0.187 L 2 = 0.492 l_ 3 = 0.316 C ] = 1.160 C 2 = 0.492 C 3 = 1.100 

Number of Iterations - 110 Computer Run-time - 670 seconds 

Termination Caused by - J^ 1+ ^ < 10 ^ 

Graphical Comparison of Desired and Actual Frequency Response 




(Hz) 



Figure 10 

Data for Example 4 
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Circuit Configuration 

h •’ M c, 

-| + 

-n. v 2 

Starting Values of the Variable Parameters 
L ] = 1.0 L 2 = 1.0 C 1 = 1.0 C 2 = 1.0 

Final Values of the Variable Parameters 
L 1 = 0.301 L 2 = 0.745 C ] = 0.877 C 2 = 0.345 

Number of iterations - 48 Computer Run-time - 187 seconds 

Termination caused by - J^ 1+ ^ < 10"^ 

Graphical Comparison of Desired and Actual Frequency Response 
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Circuit Configuration 




Starting Values of the Variable Parameters 

L, = 2. L 0 = 2. L, = 1. L = 1 . C, = 2. C„ = 2. C 0 = 1 . C =1. 
1 2 3 5 1 2 3 4 

c 5 -i. 

Final Values of the Variable Parameters 
L 1 = 1.551 L 2 = 1.762 L 3 = 1.452 L 4 = 1.043 L g = 0.697 

C ] = 1.889 C 2 = 1.750 C 3 = 1.355 C 4 = 0.948 C 5 = 0.374 

Number of Iterations — 99 Computer Run-time — 1930 seconds 

Termination caused by - J^ 1+ ^ < 10”^ 

Graphical Comparison of Desired and Actual Frequency Response 




Desired — 
Actual 



f 

(Hz) 
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run-time of nineteen minutes). This demonstrates the price paid to 
achieve a higher degree of "closeness of fit" of the actual to the 
desired response, 

C. PROBLEM AREA 

In attempting to determine the maximum capability of the program, 
a complex sixteen-element low-pass filter was attempted. Results showed 
an unusually long computation time for one iteration (approximately 200 
seconds per iteration). Investigation disclosed that about eighty per- 
cent of the computational time was spent in the tree-finding subroutine 
of the modified Calahan section. As a result, after using 720 seconds 
of computer run-time and obtaining only 3.25 iterations, the problem was 
abandoned as requiring excessive computation time for a reasonable solu- 
tion. 
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V. A GUIDE TO THE USE OF THE OPTIMIZATION PROGRAM 



In this section the general characteristics of the optimization 
program are first discussed, followed by examples which are explained in 
detail and illustrate the use of the program. Finally, several guidelines 
are offered that result from experience with the program. 

A. GENERAL CHARACTERISTICS 

The optimization program consists of three main parts: 

1. The gradient-projection search section 

2. The modified Calahan program, which produces the frequency response 
of the given network. 

3. Subroutine GRAD which acts as the interface between sections one 
and two, and computes values for the gradient vector and the cost 
function. 

The program is dimensioned to accept fifteen variable network para- 
meters in a circuit with as many as 100 passive elements and 30 nodes. 

The circuit to be optimized is topologically described using nodes, 
branch element types, and their values. For the variable parameters, the 
specified value is the "starting point"; for elements that remain fixed 
the actual value is specified. Provisions are made for use of constraints 
and a cost-function weighting matrix, h. Termination criteria and the 
desired frequency response are supplied by the user. Program output 
includes: final parameter values; final values of the cost function and 

the norm of the gradient; the total number of iterations; printer plots 
of frequency response, phase, and delay for the final element values. 

Total storage space required for the computer program is 120K bytes. 

B. EXAMPLES 

Figure 13 is a general explanation of the required and the optional 
data cards necessary for use of the program. The specific format for each 
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1 2 3 4 5 

X I / - /- / - 



6 7 




I la 

I I b 

III 



IV 



V 

VI 

VII 

VIII 



1. NVAR — the number of elements to be varied 

2. NCON — if NCON is set to 0, Card I la is 
used to designate common upper and lower 
bounds; if NCON is set to the number of 
constraint equations. Card group lib is 
used, and Card I la is omitted. 

3. KWIT - Maximum number of iterations desired 

4. EPS4 — e of equation (2.3) 

5. WM(1) — h^ in the weighting matrix, and 
indicates that Card group VIII will be 
provided; if hi is set equal to 0, equal 
weighting of all frequency response points 
will be used by the program 

6. The common lower constraint 

7. The common upper constraint 

8. This block of cards designates all con- 
stants in the constraint matrix, if required 

9. The number of passive elements 

10. The number of active elements 

11. The number of nodes 

12. The positive input node number 

13. The negative input node number 

14. The positive output node number 

15. The negative output node number 

This block of cards describes the circuit 

topology. The first n cards must be the 

elements to be varied. 

16. The "from" node number 

17. The "to" node number 

18. The element type (R, L, C) 

19. The element value 

20. The number of frequency intervals 

21. The minimum frequency (Hz) 

22. The maximum frequency (Hz) 

23. This block of cards describes the desired 
frequency response in decibels 

24. This block of cards is required only if 
the user desires to specify the weighting 
matrix, and is necessary if item (5.), is 
greater than 0 



Figure 13 

Data Cards Necessary for Use of the Optimization Program 
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card (or group of cards), is shown in the two examples that follow. It 
should be noted that two of the entries in Card I serve dual functions. 

1. If the second entry (NCON, which indicates the number of con- 
straint equations) is set equal to zero, upper and lower bounds for all 
variable elements are specified with Card I la , satisfying the equation 

a <_ p^ <_ b , i = l,...,n where a is the common lower bound and b is the 
common upper bound; if a positive integer is entered. Card group lib 
must be provided by the user, specifying the constants in the constraint 
equations, omitting Card I la. (The formation of the constraint matrix 
is discussed in subparagraph B.2. of this chapter.) 

2. If the fifth entry in Card I is set equal to zero, all frequency- 
response points will be weighted equally, and Card group VIII is omitted; 
if any positive number other than zero, is entered, the program stores this 
value as the first element of the weighting matrix, h-j. Card group VIII 
(containing K-l cards) representing h 2 »...,h^ must then be provided as 
the last section of the data deck. 

Two examples will be given using a four-element filter. The first 
example demonstrates use of the program in the simplest fashion, while the 
second illustrates the use of more complex constraint and weighting matri- 
ces. 

1 . Example 1 

The circuit to be optimized is shown in Figure 14(a). The load 
resistor, R, is to be held constant; all other elements are to be opti- 
mized to match a given frequency response. A maximum of fifty iterations 
and e = 10" 6 (see equation (2.3)) are the termination criteria. All 
frequency points are to be equally weighted; eleven frequency points will 
be matched over the frequency interval (0.10,0.20). The three variable 
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CARD NUMBER 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 



DATA 

3 0 50 .000001 0. 

0.1 4. 

4 0 4 1 4 3 4 

1 2 L 2.0 

2 4 C 1.5 

2 3 L 1.5 

3 4 R 1.5 
10 

0.10 0.20 

-0.005018 



FORMAT 

(3I4.2E15.8) 

(2F10.0) 

(7 Cl 2 , 1 X ) ) 



► (2(I2,1X),A1 ,1X,F10.0) 
(13) 

(2E15.8) 

► (E15.7) 



20 



-10.78127 



(b) 

Figure 14 



Data Cards for Example 1 
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parameters are to be constrained by the inequality 0.1 <_p^ <_4.0, 
i = 1,2,3. Figure 14(b) lists the required data cards. An explanation 
of the information on the cards follows. 

Card 1 . In sequence: 

3 - The number of variables, n. 

0 - The number of constraints; in this case, 0 signifies 
that all variables are limited by like constraints, 
the values of which are specified on the next card. 

50 - The maximum number of iterations. The program will 
terminate with normal output after the 50 Ln iteration, 
if no other termination criteria has been reached. 
0.000001 - This designates e. 

0. - Assignment of zero to this field indicates that all 
frequency-response points are to be weighted equally; 
that is, h. = 1., i = 1,...,K. 

Card 2. This designates the lower constraint as 0.1 and the upper 
constraint as 4.0 for all variables. 

Card 3. In sequence: 

4 - The total number of passive circuit elements 

0 - The total number of active circuit elements 
4 - The number of nodes 

1 - The positive input node number 
4 - The negative input node number 

3 - The positive output node number 

4 - The negative output node number 

Cards 4-7. These cards describe each branch element, as specified 
earlier. Note that the card for the non-variable load 
resistor follows the cards for the three variable elements. 
Card 8. The number of frequency intervals (the number of frequency 
points to be matched minus one) is entered. 

Card 9. The lower and upper frequency limits are entered. 

Cards 10-20. These cards list the desired frequency-response values 
at equal intervals, from the lowest frequency to the highest. 
There must be K points specified. The values must be in 
decibels. A maximum of 100 frequency points may be entered. 
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CARD NUMBER 



DATA 



FORMAT 



1 

2 

3 

4 

5 



3 7 100 „ 000001 

1. 

0 . 

Oo 

1 . 



22 


1. 




23 

0 


J 

0 




0 

29 


0 

0 

Oo 




30 


4 


0 4 1 


31 


1 


2 L 2o0 


32 


2 


4 C 1.5 


33 


2 


3 L 1 o 5 


34 


3 


4 R 1.0 


35 


10 




36 


.1 




37 


-o005018 

0 


47 


0 

-10o78127 


48 


2. 




49 

0 


3. 

0 




0 

57 


0 

lo 





1. (3I4,2E1 5„8) 

>(E18.8) 

(7(12, IX)) 

•v 

>(2(I2,1X),A1,1X,F10.0) 

(13) 

(2E15o8) 

>(E15o7) 

> (El 5,8) 



Figure 15 

Data Cards for Example 2 
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2. Example 2 



The same circuit is to be optimized, in this case using specific 
constraint and weighting matrices. The constraints are: 

0.1 <_ p^ £ 1.5 
0.5 <_ $2 £ 2.0 
1 .0 £ p 3 £ 3.0 
2p 2 < P 3 

These equations must be put into the form 

Ap - B £ 0 * (5.1) 



Proper manipulation of the equations yields the seven constraint 
inequalities shown below in matrix form. 



1. 


0 . 


1 

0 . 




0.1 


0 . 


1. 


0 . 




0.5 


0 . 


0 . 


1. 




1.0 


-1. 


0 . 


0 . 


p - 


-1.5 


0 . 


-1. 


0 . 




-2.0 


0 . 


0 . 


-1. 




-3.0 


0 . 


-2. 


1. 




0 . 



The weighting matrix to be used is 



h 



[1. 2. 3. 4. 5. 6. 5. 4. 



3. 2. 1.] 



This weighting puts greater emphasis on a close fit at the center portion 
of the frequency response curve. The number of elements in the weighting 
matrix must correspond to the number of frequency points, K. Figure 15 
lists the required data cards. An explanation of the information on the 
cards follows. 



This notation means that each component of the vector Ap - B is 
non-negative. 



Card 1. In sequence: 



3 - The number of variables , n« 

7 - The number of constraints. This corresponds to 
the number of rows in A, and cannot exceed 32 
without re-dimensioning the program. 

100 - The maximum number of iterations. 

0.000001 - This designates e. 

1. - This assigns a positive value other than zero to 
the first element of the weighting matrix, and 
indicates that specified weighting matrix values 
will be provided. 



Cards 2-22. These cards list the values of each element of 
the constraint matrix A, sequentially by row. 

Cards 23-29. These cards list each value of the constraint 
matrix B. 



Cards 30-47. No change from example 1. 

Cards 48-57, These cards list the ten remaining values of 
the elements of h. 



C. GUIDELINES 



Experience permits the suggestion of several guidelines that may be 
helpful in using the program. 

1 . Termination Criteria 

Appropriate values for the maximum number of iterations and for e, 
depend primarily on the closeness of fit desired, Higher values for the 
maximum number of iterations (100 or higher), coupled with lower values 
of e (10 or smaller) will ensure a closer fit, but at the expense of 
computer run-time. 

Additionally, the value for Y in equation (2.4) 

ao(p) 



< Y 



ap 
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is designated within the main program as "ECON" and is set at a value of 
-4 

10 . This termination criterion may be altered by the user if desired. 

2. Closeness of Fit at Termination 

Closeness of fit of the final frequency response to the desired 
response is somewhat dependent on how closely the circuit is capable of 
approximating the desired frequency response. If a close enough fit is 
not achieved, the following steps are recommended. 

Step 1. Impose more stringent termination criteria and re- 
initialize with the final values just obtained. (This step should be 
omitted if it is clear that the circuit configuration is not capable of 
coming close to the desired response.) 

Step 2. Try different starting points for the variables if the 
termination value of the cost function is high. 

Step 3. Try a different network configuration. 

In addition, the value N in equation (2.2), the cost-function 
equation, has been set at N = 2 for least-squares curve fitting. This 
may be changed to N = 4 to penalize more heavily large deviations from 
the desired response. Statements 99 and 101 in subroutine GRAD both 
contain "**2" which should be changed to "**4" to accomplish this alter- 
ation. 

3. Number of Frequency Points 

Enough frequency points should be specified to completely describe 
the desired response. Doubling the number of frequency-response points 
does not double the computer run-time, but requires an increase of only 
about five to twenty percent in time (the larger percentage applies to 
small circuits, the smaller percentage to large circuits). This is 
because most of the run-time used is spent producing the network function 
(poles and zeros); once this is done, computation of the frequency 
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response is relatively fast. Failure to adequately describe the desired 
response may produce poor results. Experience with a limited number of 
examples indicates that an adequate description is provided when fre- 
quency points are provided to at least ten decibels below the desired 
flat response for low- and band-pass filters, and spaced over a close 
enough interval so that "corners" are well-defined, 

4, Computation Time 

Complex circuits (those with many trees) and any circuit with 
more than about eight elements require run-time in excess of ten minutes 
(on the IBM-360 digital computer). Investigation has shown that this 
time requirement is a function that varies directly with the number of 
trees in the circuit. Appropriate alteration of the Calahan subroutines 
to improve or alter the tree-finding method should provide a decrease in 
computation time necessary to produce a reasonable solution for large 
circui ts. 
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VI. COMPARATIVE RESULTS OF GRADIENT-PROJECTION AND 
PATTERN-SEARCH METHODS 



Optimization of filter frequency response using pattern-search 
methods is the subject of a Naval Postgraduate School thesis being 
concurrently written by Captain James Lau, USMC. A comparison of the 
results of attempts to optimize the frequency response of three networks 
using each of the programs is shown in Figure 16. The first and second 
comparisons were optimized to a totally realizeable response, while the 
third was to a straight-line approximation. In all cases, the pattern- 
search optimization program proved faster and capable of producing more 
accurate results. 
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GRADIENT PROJECTION 


PATTERN SEARCH 




FILTER TYPE 




RESULTS 


RESULTS 


1 . 


Eight-element band-pass 










* 

Final cost-function 


value 


Oo 564 x 10“ 4 


0.220 x 10" 5 




Computation time 




328 seconds 


184 seconds 


2. 


Five-element low-pass 










Final cost-function 


value 


0.458 x 10" 2 


0.204 x 10" 5 




Computation time 




133 seconds 


49 seconds 


3. 


Five-element low-pass 










Final cost-function 


value 


1.655 


1 .651 




Computation time 




132 seconds 


50 seconds 



Figure 16 

Results of Three Network Optimizations, 

Comparing Gradient-Projection and Pattern-Search Methods 
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VII. CONCLUSIONS 



The purpose of this investigation was to determine the feasibility 
of using minimization techniques to cause a given network configuration 
to match a desired frequency response as closely as possible. This 
method was to use the computer's speed coupled to optimization techniques 
to provide the filter designer with a tool that quickly determines the 
network element values that most accurately cause the desired frequency 
response to be reproduced. Two existing programs capable of 

1. minimizing a given function, and 

2. producing a given network frequency response 

were adapted to produce a program capable of accomplishing the above- 
stated purpose. 

The results presented in Chapter IV and the comparisons in Chapter 
VI lead to the following conclusions: 

1. It is feasible to obtain realistic approximations to a 
given frequency response, using a circuit capable of attaining that 
response. 

2. Satisfactory solutions using gradient-projection techniques 
can be obtained for filters with fewer than about twelve elements. An 
unreasonable amount of time may be required in attempting to optimize 
larger circuits. 

3. The pattern-search method for minimization is shown to be 
superior to the gradient projection method in both time and accuracy, for the 
examples compared. 
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Recommendations for logical extensions of this investigation follow: 
1„ Investigate a means to determine the gradient vector more 
accurately and more quickly. 

2. Alter the analysis program so that the tree-generating 
procedure is used only when necessary, rather than in each iteration. 
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APPENDIX A 

PROCEDURE FOR FINDING THE PARABOLIC MINIMUM USED IN THE 
FUNCTION-SEARCH INTERPOLATION METHOD 



I. The general equation of a parabola is: 
Ax 2 + Bx + C = y 



(A.l) 



The unknown coefficients A, B and C that will cause a parabola to pass 
through any three known points (not in a straight line) can be found 
by solving three simultaneous equations: 



Ax^ 2 + Bx^ + C = y 

Ax^ 2 + Bx^ + C = y 2 
2 

AXg + BXg + C = y^ 

2. These equations can be solved using the matrix equations: 




(A. 2a) 
(A. 2b) 
(A. 2c) 



(A. 3) 



3. The minimum of this parabola can be found by setting the first 
derivative of (A«l) with respect to x equal to zero and solving for x . . 



3x ( Ax 2 + 6X + C - y ) - 2Ax m . n ♦ 8 ■ 0 



X 



min 



_B_ 

2A 
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APPENDIX B 



FUNCTION-SEARCH INTERPOLATION EXAMPLES 



1. Shown below are two examples (many others were noted, but are not 
included herein) where the gradient-projection minimization process 
seemed to be approaching a minimum, (no conclusive proof of this is 
offered, only the indications listed herein) when a distinct change in 
all of the variables as well as the cost function occurred. In the first 
four of the five iterations shown, the changes in the parameter values 
were small, and the cost function was being very slowly reduced. This 
behavior is characteristic of gradient methods when a minimum has been 
found (this can also be representative of moving down a "ridge line" or 
a "narrow valley"). This "jump" would not take place if the gradient 
method (as opposed to the gradient-projection method) were used. (Only 
five iterations are shown, as those before and after are not pertinent 
to the illustration.) 

First Example. Five-element (four-variable) low-pass filter 



Iteration 


Variable Parameter 


Values 


Cost 


number 


Pi 


P 2 


P3 


P 4 


function 


19 


1 .3311 


1 .9541 


0.9654 


0.5633 


0.2425 


20 


1 .3284 


1.9489 


0.9624 


0.5627 


0.2363 


21 


1.3523 


1 .9178 


0.9801 


0.5493 


0.2171 


22 


1.3468 


1.9108 


0.9743 


0.5489 


0.1978 


23 


1 .4730 


1 .6550 


1 .0565 


0.4367 


.02567 



Second Example. Six-element (five-variable) low-pass filter 



Iteration 




Variable 


i Parameter Values 


Cost 


number 


p l 


p 2 


p 3 


P 4 


P5 


function 


32 


1 .7525 


1 .8321 


1 .4150 


1.0029 


0.4440 


0.02815 


33 


1.7493 


1.8270 


1 .4163 


1 .0041 


0.4469 


0.02782 


34 


1.7503 


1 .8269 


1.4185 


1 .0037 


0.4488 


0.02733 


35 


1.7497 


1 .8263 


1.4182 


1.0039 


0.4503 


0.02730 


36 


1 .6674 


1.6507 


1 .5231 


1.0256 


0.5055 


0.01116 



54 



2. These examples tend to substantiate the conclusion that the function- 
search interpolation method used in the gradient projection minimization 
program increases the likelihood of converging to the absolute minimum, 
for when the vector in the negative gradient direction intersects a 
point with a lower cost-function value — perhaps relatively distant from 
the base point — this new smaller-valued location will be taken up as 
the new base point. 
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6 1 DROP, IRAX,NAC, IDEP.NACMl , JV,NMACT ,NQP1 , NS , NS PI , NVAR 2 , INC 
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GRID IS THE MATRIX USED TO PLOT THE POINTS 
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